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OPTIMAL CURE CYCLE DESIGN 
OF A RESIN-FIBER COMPOSITE LAMINATE 

Jean W. Hou and Jeenson Sheen 
Old Dominion University, Norfolk, VA 23508 

ABSTRACT 

High performance polymeric composites have been experiencing 
increasing usage in the aerospace and automobile industries. Such 
materials are commonly composed of long or chopped fibers embedded in 
the thermosetting resin matrix. Changes in physical and chemical 
properties of such composite materials during the cure process are 
rather complex. Thus, it is not a trivial task to properly design a 
cure cycle (temperature and pressure profiles) for a cure process. The 
material should be cured uniformly and completely with the lowest void 
content; the temperature inside the laminate must not exceed some 
maximum value; and the cure process should be completed within the 
shortest amount of time. In the past, most cure cycle designs for newly 
developed composite systems are based upon the technique of trial and 
error. Such approach has long been recognized as costly and 
inefficient. Several simulation models have been developed recently for 
curing various epoxy matrix composites. This development represents a 
significant advancement in computerizing the cure cycle design. The 
next quest comes naturally to be the search of the "best" cure cycle for 


* Assistant Professor, Department of Mechanical Engineering and 
Mechanics 

** Graduate Student 


i 


a given composite laminate. The major thrust of this thesis is to study 
a unified Computer-Aided Design method for the cure cycle design that 
incorporates an optimal design technique with the analytical model of a 
composite cure process. The preliminary results of using this proposed 
method for optimal cure cycle design are reported and discussed in this 
thesis. 

The cure process of interest is the compression molding of a 
polyester which is described by a diffusion-reaction system. The finite 
element method is employed to convert the initial-boundary value problem 
into a set of first order differential equations which are solved 
simultaneously by the DE program. The equations for thermal design 
sensitivities are derived by using the direct differentiation method and 
are solved by the DE program. Finally, a recursive quadratic 
programming algorithm with an active set strategy called a linearization 
method is used to optimally design the cure cycle, subjected to the 
given design performance requirements. The difficulty of casting the 
cure cycle design process into a proper mathematical form is recognized 
in this study. Various optimal design problems are formulated to 

address these aspects. The optimal solutions of these formulations are 
compared and discussed, and the major parameters which play major roles 
in the cure cycle design for a given composite laminate are identified. 
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CHAPTER 1 
INTRODUCTION 

Continuous fibers reinforced polymeric composites are used in many 
applications where high strength and low weight are important, such as 
in the aerospace industry. These composites are produced by embedding 
long or chopped high strength reinforcing fibers in a polymeric 
matrix. The matrix is important not only for transferring the load 
between fibers, but also for providing resistance to fracture in the 
composite and to distortions caused by the environment. One means to 
manufacture the composite laminate with continuous fiber is by combining 
unidirectional ly oriented layers of fibers pre-impregnated with the 
uncured resin matrix and subsequently laying these up in the desired 
directions. This pre-impregnated composite precursor is called a 
prepreg. The composite parts and structures are then manufactured by 
curing the prepreg material. The cure process is accomplished by 
exposing the prepreg material to the elevated temperatures and pressures 
for a predetermined length of time. These elevated temperatures and 
pressures to which the prepreg material is subjected are referred to as 
the cure temperature and the cure pressure. The cure temperature 
provides the heat required for initiating and maintaining the chemical 
reactions in the resin during the cure process. Once the reactive 
matrix resin is heated, polymerization reactions start to occur. 
Polymer chains are extended and later crosslinked. Interactions between 
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the heat generated from the reactions and the heat transmitted by 
conduction and convection from the surrounding environment create a 
highly non-linear chemoviscosity profile which can dictate the resin 
flow characteristics during the cure process. After the resin is 
melted, by applying the cure pressure, the excess resin is squeezed out 
of the material. Consequently, the individual plies are consolidated 
and the vapor bubbles are compressed. 

The magnitudes and durations of the cure temperature and the cure 
pressure applied during the cure process have significant consequences 
on the performance of the finished product. It is not a trivial task to 
design a proper cure cycle, that is, to determine a cure temperature and 
a cure pressure, due to the complex nature of the changes in physical 
and chemical properties of the composite material during the cure 
process. The composite material should be cured uniformly and 

completely with the lowest void content; temperatures inside the 
laminate must not exceed some maximum value; and the curing process 
should be completed in the shortest amount of time [1]. In the past, 
the cure cycle is selected empirically for newly developed composite 
systems. This trial and error type of empirical approach has long been 
recognized as costly and inefficient. The shortcomings of the empirical 
approach could be overcome, however, by the use of an analytical 
model. The recently developed model for curing epoxy matrix composites 
by Loos and Springer [1] represented a significant advancement in this 
aspect. The cure cycle can now be designed in a systematic manner to 
meet the major considerations mentioned earlier for various composite 
resin/fiber systems. As for the composite with chopped fibers, its 
manufacturing process [2] may be different from the one stated above. 
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However, the process can still be characterized similarly by the resin 
flow, chemical kinetics, chemoviscosi ty, etc. The considerations in the 
selection of the proper cure cycle are also very similar. 

1.1 Objective and Motivation 

As mentioned above, currently the method for the cure cycle design 
is a type of parametric study conducted by either physical experiments 
or numerical simulation [1,2], This type of approach is costly and 
inefficient; and the true optimal condition may never be found. 
Therefore, the development of a unified and systematic approach to 
obtain an "optimal" cure cycle for a given composite material is 
needed. In this study, a unified Computer-Aided Design (CAD) method is 
proposed, which incorporates an optimal design technique with an 
analytical model of a composite cure process. 

The optimization technique is not new and has been successfully 
applied to many transient problems in various disciplines such as 
mechanical systems [3], structural dynamics [4], chemical process 
control [5], and optimal control [6]. However, to the author's 
knowledge, the application of optimization technique to the cure cycle 
design has never been studied before. The goal of this study is to 
exploit the potential of applying the optimization technique to the cure 
cycle design. 

As an initial attempt at this goal, this study is concentrated on 
the optimal design of the cure temperature only. The press molding 
process of a polyester [2] is used as an example. Various optimization 
formulations for the cure cycle design are suggested. Their optimal 
solutions are obtained by using the integrated computer program which is 
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developed for the simulation of a heat conduction model coupled with a 
chemo-kinetic model during the process and for the optimization of the 
cure cycle design. Those optimization formulations are evaluated based 
on their performances and merits. The experience gained from this study 
is very important and valuable for the further application of the 
proposed CAD approach to a more realistic composite manufacturing 
process. 

1.2 Scope of the Present Work 

The major contents of this thesis are Thermal Analysis, Thermal 
Design Sensitivity Analysis, and Optimal Cure Cycle Design. The 
considerations of the computational formulation of the reaction- 
diffusion system are reported in the next chapter. The finite element 
discretization is introduced to convert the initial-boundary value 
problem into a set of first order differential equations. These 
equations are then solved simultaneously by the DE program [7]. To 
verify the computational procedure, the numerical results of examples 
are compared with the analytical solutions or the solutions found in the 
literature [2,8,9,10]. 

In Chapter 3, two methods, the adjoint variable method and the 
direct differentiation method, have been studied for the thermal design 
sensitivity analysis. It is concluded that the direct differentiation 
method is superior to the adjoint variable method in terms of accuracy 
and physical interpretation. Once the information of design sensitivity 
is provided accurately, any gradient-based mathematical programming 
algorithm can be applied to generate the optimal design iteratively. 
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In Chapter 4, an iterative scheme is outlined to determine the best 
cure process. Numerical calculations are carried out by using the 
finite element method for analysis and a recursive quadratic programming 
technique for optimization. The results of the optimal cure cycle 
design for various formulations of objective functions are also 
presented. 

Finally, remarks and conclusions of the present work are given in 
Chapter 5. 
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CHAPTER II 
THERMAL ANALYSIS 

One major step of manufacturing composite materials which consist 
of thermosetting resin matrix is using the elevated temperature and 
pressure to cure the materials. The analytical model of such a cure 
process is rather complicated. Such an analytical model should be able 
to address: 

a) the heat transmission between the surrounding environment and 
the material, 

b) the temperature distribution within the material, 

c) the heat generated by the chemical reaction of the resin, 

d) the phase change of the resin characterized by its 
chemoviscosi ty, and 

e) the non-newtonian flow problem associated with the cure pressure 
and the liquidized resin. 

Among others, the major parameters which play important roles in 
modelling the cure process may be identified as the temperature, the 
degree of cure of the resin, the resin viscosity, the resin flow 
velocity, the thickness of the processing material and the material 
properties, such as density, heat conductivity, etc. Most parameters 
are varied with respect to the time and spatial location. 

To focus on the optimization aspects of study, several assumptions 
which simplify the formulation of the cure process are introduced. One 
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is that the temperature and the degree of cure of the composite material 
are uniformly distributed in a plane parallel to the tool plates. As a 
result, the parameters of the problem are functions of the time and the 
position along the thickness of the composite only. Secondly, the resin 
content in the material is assumed to be low so that only a small amount 
of resin will be squeezed out of the material. Consequently, the resin 
flow model as well as the pressure cycle can be neglected from 
considerations. Furthermore, the material properties which depend on 
the fiber-resin ratios can be assumed to be constant as well. Finally, 
it is assumed that there is no deviation between the cure temperature 
and the temperature on the surfaces of the material adjacent to the tool 
plates. Thus, the cure temperature can be considered as the boundary 
temperature of the material. The temperature distribution and the 
degree of cure are the only two state variables which are governed by a 
diffusion-reaction system. After simplification, the diffusion-reaction 
system can be used to describe the cure process of a continuous fiber 
composite with low resin content [8], and the cure process of a 
composite with the chopped fibers [2]. The interest of this study, 

however, is in compression molding of a filled polyester resin 
reinforced by chopped glass fibers. The unmolded composite is produced 
in sheets which are from 3 to 6 mm thick, typically. The resin consists 
of a thick dough and the chopped fibers (about 25 mm long) which are 
randomly oriented in the plane of the sheet. In this form, the material 
is called sheet molding compound, or simply SMC. 
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2.1 Diffusion-Reaction System 

A diffusion-reaction system consists of an one-dimensional unsteady 
heat-conduction equation with its boundary and initial conditions, and a 
kinetic model of the cure rate of the resin. By using the finite 
element method and a numerical integrator, a computational algorithm is 
developed to analyze the thermo-chemical reaction. The numerical 

results of several examples are obtained by the developed scheme and 
compared with the existing solutions. 

The temperature distribution T(z,t) and the degree of cure a(z,t) 
of the resin inside the composite depend on the rate which heat is 
transmitted from the environment into the material. The temperature 
inside the composite can be calculated by using the law of conservation 

of energy together with an appropriate expression for the cure 

kinetics. By neglecting energy transferred by convection, the energy 
equation may be expressed as 

pcT = kT" + pHra (2.1) 

with the boundary conditions, 

r ( 0, t) = 0, 0<t<x, (2. 2. a) 

T(h,t) = T c (t) 0<t<-c, ( 2 . 2 . b ) 

and the initial condition, 

T(z,0) = T q (z), 0<z<h (2.3) 

where p and c are the density and the specific heat of the composite 
material, respectively, k is the thermal conductivity in the direction 
perpendicular to the plane of composite material, and 2h is the total 
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thickness. The temperature and the total or ultimate heat of reaction 
during the cure process are denoted by T and Hr, respectively. 
Moreover, the dot on the top of a symbol indicates the time 
derivative and the superscript denotes the spatial derivative. 
According to the assumptions discussed previously, all coefficients in 
Eqs. (2.1) - (2.3) are treated as constants. Note that the cure 

temperature T c (t) appears in the equation of the boundary condition Eqs. 
(2.2); and the last term in Eq. (2.1), pHra, is the rate of heat 
generated by the chemical reaction characterized by the state of cure 
a. 

The state of cure depends strongly on the resin temperature. An 
empirical expression is usually employed to address the relation between 
the state of cure and the resin temperature. The degree of cure is 
defined as 

a = H ( t ) /Hr (2.4) 

where H(t) is the heat evolved from the beginning of the chemical 
reaction of resin to some intermediate time t. Both H(t) and Hr in Eq. 
(2.4) can be measured experimentally by, as an example; Differential 
Scanning Calorimetry (DSC). For an uncured material, a approaches to 
zero, and for a completely cured material, a approaches to one. For 
example, the cure rate equation of a stepwise isothermal cure process, 
can be defined for a polyester as follows 

a = f(a,T) 

= (Kj+ K 2 a m )(l - <x) n 
= (a 1 e" d l/RT + a 2 e' d 2/RT a m )(l - a) n 


(2.5) 
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where a^, d^, d^, m, and n are constants, R is the universal gas 

constant, and and l <2 are exponential functions of the temperature. 
This example is taken from the research of compression molding conducted 
at the General Motors Research Laboratory [2], 

Some observations of interest are stated as follows: 

1. The state equations of the cure process are coupled with two state 
variables, namely, the temperature distribution T(z,t) and the degree of 
cure a(z,t) . 

2. The non-homogeneous boundary value, T c (t), is to be considered as 
the design variable. 

By using the following replacement of the temperature T(z,t) 
as 

T(z,t) = T(z,t) + T c (t), (2.6) 

the heat-conduction problem stated in Eqs. (2.1) - (2.3) can be 

simplified as an equation of T(z,t) 

pcT = kT - pcT + pHrf (a,T ,T ), (2. 7. a) 

C 0 

with the homogeneous boundary conditions, 

a 

T (0,t) = 0, 0<t<T, (2.7.b) 

T(h,t) = 0, 0<t<x, (2.7.C) 

and the initial condition, 

T(z,0) = T (z) - T (0), 0<z<h (2.7.d) 

0 c 

where f is defined in Eq. (2.5). From here on T is abbreviated as T for 
further simplification. Because the initial temperature T 0 (z) is the 
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same as the initial cure temperature for most applications, the Eqs. 
(2.7) may have a homogeneous initial condition as well. Note also that 
the cure kinetics of the polyester are expressed in terms of the 
absolute temperature. Thus, the Kelvin's degree is used in this study 
as the unit of temperature. 

2.2 Finite Element Model 

The finite element discretization is introduced herein to convert 
the initial-boundary value problem, Eqs. (2.5) and (2.7), into a set of 
first order differential equations. These equations are then solved 
simultaneously by a numerical integration code called DE [7]. 

The quadratic and linear polynomials are used to interpolate the 
states of the temperature and the degree of cure, respectively. The 
notations ^i ’ ancl a i denote tfie temperature, the temperature 

gradient and the degree of cure, respectively, at node i as shown in 
Fig. 2.1. 

There are two degrees of freedom assigned at each node to 
approximate the temperature distribution. The interpolation functions 
of the temperature in each element are then defined as follows 

N 2 (z) =1-3 z 2 /* 2 + 2 z 3 /! 3 

N 2 (z) = -z + 2 z 2 /l - z 3 /A 2 

( 2 . 8 ) 

N 3 ( z ) = 3 z 2 /l 2 - 2 z 3 /A 3 

N 4 (z) = z 2 /i - z 3 /A 2 


Hence the 
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where the notation of A is the length of an element, 
temperature within an element is then given as 


T(z,t) = {Nj(z) , N 2 (z), 


N 3 (z), N 4 (z)} { 


Tj ( t) 

T 2 (t) 
T 3 (t> 
T 4 (t) 



(2.9) 


As for the degree of cure, each node has one degree of freedom. The 
interpolation functions are defined as 


!_! (z) = 1 - z/A 
L 2 (z) = z/A 


( 2 . 10 ) 


where A is the length of an element. Hence, the degree of cure within 
an element is defined as 


a(z,t) = {LjU), L 2 ( z) } 


al ' ' 

( 1 

a 2 ( t) 


= L 



( 2 . 11 ) 


Note that the shape functions employed here for the temperature 
distribution and the state of cure are the same as those for an elastic 
beam and a truss, respectively. 

Equations (2. 7. a) and (2.5) provide integral identities for any 
arbitrary testing functions u(z) and v(z): 


0 = J (pcT - kT + pci - pHrf) udz 

c 

o 


= / (pcTu + kT u + pcT u - pHrfu) dz - kT u 

c o 

o 


NE i+1 

= E / (pcTu + kT'u' + pcT u - pHrfu) dz 
1 z. C 

i 


( 2 . 12 ) 


1 1 

= / (a - f) v dz 


NE Z i+1 


- f) v dz 


(2.13) 


Where the testing function u(z) has to satisfy the kinematic boundary 
condition, Eq. (2.7.c) and NE denotes the number of elements. The 
boundary terms of Eq. (2.12) drop because of u(h) = 0, and T'(0,t) = 0. 
The last equality in Eq. (2.13) indicates that the domain (0<z<h) has 
been divided into NE finite elements. Considering any interpolation 
functions Nj and Lj defined previously as the testing functions u(z), 
and v(z), respectively, and using Eqs. (2.9) and (2.11), one can 
establish two matrix equations for the nodal temperature and the nodal 
degree of cure as 


pc[C]T + k[K]T + pcPT = Q(a,T,T ), 

- c c 


(2. 14. a) 


= R(a,T,T 


(2.14.b) 


with the boundary conditions. 
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I = {0, T 2 (o), T 3 (o), •••, T 2n _ 1 (o), 0} T 


(2.14.C) 


where the unknowns are vectors T = {T^, T 2 , •••» T^}^ and a = {a^.o^, 
.... a for n = NE + 1 as the total number of nodes. Note that T,(t) 
and T 2n (t) are always zeros in the analysis according to the boundary 
conditions, T(h,t) = 0 and T' (0,t) = 0, respectively. The initial 

conditions for each of the components in T and a are 


T 


2i-l 


(0) = T 


Room 


- T c (0), 


T 21 (0) = 0, 
and 

a.(0) = 0 

where i is the nodal number, and T 2l -_ 1 and T 2l - are corresponding to the 
temperature and the temperature gradient at the node i. Moreover, T Room 
represents the temperature of the composite at the onset of the cure 
process that is usually uniformly distributed. The detailed derivation 
of the matrix equations, Eq. (2.14), is given in Appendix A. The above 
equations, coupled with each other through the non-linear terms on the 
right side of Eq. (2.14), can be solved simultaneously by using a 
computer code called DE. 

The DE program is one of predictor-corrector integration algorithms 
using Adams family of formulas. The truncation error is controlled by 
varying both the step size and the order of the polynomial approxi- 
mation. The truncation error of the solution Z n+1 at time step t n+1 is 
required by the DE program to satisfy the following relation: 


trunc < ABSERR + RELERR • Z 

n 
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where Z n is the solution of the differential equation at t n , and the 
values ABSERR and RELERR are supplied by the user. 

The OE program is quite easy to be used and has the capability to 
manage moderate stiff equations which happen commonly in the problems of 
chemical kinetics. To maintain a unified accuracy in the analysis, the 
computation of two state variables, namely, the temperature and the 
degree of cure, are subjected to the same error tolerance in this study. 


2.3 Numerical Examples 


In this section, three examples have been studied. The numerical 
results are in a very good agreement with the theoretical solution in 
example one, and with the numerical solution reported in the literature 
[2,8,9,10] in example two and three. 

Example 1: one-dimensional heat conduction with prescribed end 

temperature. 

This problem can be described by the following equation, 

pcT = kT~ (2.15) 


with the boundary conditions, 

T" (0, t) = 0, 0<t<x, 

T (h,t) = 0, 0<t<x, 


and the initial condition, 

T(z,0) = T q , 0<z<h. 

The exact solution of this problem [9] is given as 


CD 

T (z,t) = l 

n-1,3,5. . 


4T 2 2 2 

o nirz (n % kt/pcX ) 

nix X 


(2.16) 


where the notation X is the total thickness 2h. 
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The finite element problem is solved with A = 10, k = 1, pc = 1, 
and T q = 100. Twenty equal -spaced elements are used. Both the finite 
element and the exact solution are shown in Fig. 2.2 for the time grids 
t = 4.0, 8.0, 12.0, 16.0, 20.0, respectively. Note that a very good 
agreement exists between the finite element solution and the exact 
solution. 

Example 2: one-dimensional heat conduction problem coupled with a 
given chemical -kinetic model [2]. 

This example simulates the compression modelling of composite 
laminates with its thickness 10 mm. The initial (material) temperature 
is 298°K and the cure temperature is taken as a constant temperature, 
i.e., T c (t) = 423°K. This problem is described by Eqs. (2.1) - (2.3). 
The coefficients of this problem are given in Table 2.1. 

Using the finite element method to solve this problem, sixteen 
elements are used. The finite element results of the temperature 

distribution and the degree of cure compared with the results of the 
General Motor Research report [2] are given in Fig. 2.3 and Fig. 2.4, 
respectively. Solid curves are from the General Motor Research report; 
the dash lines are calculated by Eq. (2.14). This example also shows a 
good agreement between the calculated data and that of the General Motor 
Research report. 

Example 3: one-dimensional heat conduction problem coupled with a 

given-kinetic model of thick-section autoclave cured composites [8,10], 

One simulates the autoclave processing of 192 ply prepregs with 32 % 
Hercules resin content. The resin flow can be neglected in this 
example, because the resin content is low. To focus on the heat 
conduction and the chemical-kinetic models, the measured temperature on 
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the surface of the composite laminate is used as the boundary 
temperature, instead of the temperature of the cure cycle. Although the 
heat flux, induced by the heat convection of autoclave air temperature, 
is neglected, the numerical result is in excellent agreement with Loos' 
data [10], as shown in Fig. 2.5. 
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CHAPTER III 

THERMAL DESIGN SENSITIVITY ANALYSIS 

The derivatives of functionals of responses with respect to the 
design variables are often referred to as design sensitivity 
derivatives. Most general optimization algorithms require such 
derivatives which can be used to approximate constraints and to choose a 
search direction to obtain a set of improved design variables. It is 
then needed to have a reliable means to calculate the design sensitivity 
derivatives. 

Two commonly used methods for design sensitivity analysis, namely, 
the adjoint variable method and the direct differentiation method, are 
studied in this chapter. The comparisons of numerical results obtained 
by these two methods demonstrate that the direct differentiation method 
is a better choice for this study. 

3.1 Techniques of Thermal Design Sensitivity Analysis 

As discussed previously, the cure uniformity and completeness for 
the cure cycle design. However, it is difficult to cast the measurement 
of the cure uniformity and completeness in a precisely mathematical 
form. To focus on the study of thermal design sensitivity analysis in 
this section, the cure uniformity is simply represented by the least- 
squared integral of the deviation between the pointwise temperature and 
the averaged temperature as 
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t h ? h ? 

<|> = / [/ T dz - (/ T dzr/h] dt. (3.1) 

o o o 

The above functional defines a global sense of temperature uniformity 
across the thickness of the composite during the cure process. It is 
expected that the temperature uniformity yields the cure uniformity. 
Other mathematical expressions of the cure uniformity are to be 
discussed in the following chapter of optimal cure cycle design. On the 
other hand, the following pointwise inequality may be used to force the 
cure reaction to be completed at the end of the cure process. 

a(z,T ) > a^, 0<z<h (3.2) 

where t is the total time interval of the cure process and a^. is a 

constant assigned to indicate the completeness of the degree of cure. 

The numerical techniques to calculate the design derivatives of the 

d4> 0 da 

functional 6 and the state variable a, — - and — , are to be addressed 
o db db 

hereafter. Note that the design variable b is a parameter associated 
with the cure temperature profile. 

In general, there are four ways to calculate the thermal design 
derivatives, namely, the finite difference method, the Greene's function 
approach, the direct differentiation method (the behavior space 
approach), and the adjoint variable method (the dummy load method). The 
last two methods are often mentioned in the literature [3,4,11,12]. 
Both methods lead to a set of linear equations that have a similar 
structure to the original system. 

The computational efforts regarding the direct differentiation 
method and the adjoint variable method depend mainly on the numbers of 
constraints and design variables of concern. The direct differentiation 
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method requires the solution of a differential equation for each design 
variable, while the adjoint variable method requires the solution of an 
adjoint equation for each constraint. Consequently, the direct 
differentiation method is more efficient to calculate the design 
derivatives than the adjoint variable method when the number of design 
variables is smaller than the number of constraints, or vice versa. 

It is also known that the direct differentiation method provides 
equations of design derivatives which have exactly the same differential 
operator as that of the original equations. The equations of design 
derivatives can, therefore, be solved simultaneously with the original 
system equations subjected to the same numerical error tolerance. 
Furthermore, without extra efforts, the results of the direct 
differentiation approach provide the histories of design derivatives of 
functionals and state variables. This information is very useful for a 
designer to reconstruct the design space. One may check this 
information to see whether a design variable of concern attributes to 
the perturbation of the functional of concern consistently over a long 
or just a limited period of time. 

3.2 The Adjoint Variable Method 

The problem of interest is to derive a set of equations for 
computing the design sensitivities of a functional and the state 
variables. The variation of a function 4> with respect to a design 
variable b is defined as 



4>(t,b+e6b) - (|j( t,b) 
de 
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64 ) 


lim 

£->■0 


d4> (t,b+e6b) 
de 


e=o 


4)" 6b 


(3.3) 


where 6b is the perturbation of design variable and <\>' is defined as 
The variations of state variables can be defined by a similar 
fashion. Using the functional defined in Eq. (3.1) as an example, the 
variation of this functional is then given as 


t h 

64) = / [/ 2T - 

0 0 


2 

T> 


h 

f 


o 


T dz] 6T dz dt. 


(3.4) 


The governing equations of the diffusion-reaction system defined in 
Eqs. (2.5) and (2.7) are 

pcT = kT" - pcT + ptirf(a,T,T ), 
c c 


a = f (a,T,T ). 

c 

For arbitrary functions \(z,t) and s(z,t) it is evident that 
t h 

0 = f f X(pcT + pcT c - kT" - pHrf) dz dt, (3.5) 


t h 

0=J J S(a - f) dz dt. 
o o 


(3.6) 
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The variations of these equations yield 
x h 


df 


0 = / / (pc6TX + pc6T X - p6T~X - pHr — 6aX 

c da 

0 0 


- pHr — 6TX - pHr — 6T X) dz dt, 
oT oT c 

c 

and 

T ^ af ?>f 

0 = f I (S6a - S — 6a - S — 6T - S 5T ) dz dt. 

oa oT oT C 

0 0 c 


Integrating by parts, the above equations can be rewritten as 


x h 


0 = / / [(- P cX-kX~-pHrX |I)6T + (-pcX-pHrX 6T 

oT oT c 

oo c 


df n * 

pHrX — 6al dz dt + / kX' 6T 
oa 


dt 



0 


0 

X 

h h 

X 

h 

j kX6r 

dt + / pcX6T 

dz 

+ / pcX6T 

0 

0 0 

0 

0 


dz. 


and 


x h 

0 = fj (-$6a - S 6a - S 6T - S 6T ) dz dt 

oa oT oT C 

0 0 C 


+ / S 6a 
o 


dz. 


(3.7) 


(3.8) 


(3.9) 


(3.10) 


Adding Eqs. (3.4), (3.9), and (3.10) up, one has the variation of the 
functional <|> as 
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t h h 

6* « J J [ (2T- | / Tdz) - pcX - kX~ - pHr |I -S |I]6T dz dt 
n o r o r 

oo o 


+ / / (- pHrX — - S - S — ) 6a dz dt 

oa 9a 

0 0 


t h 

+ J J ( -pcX - pHrX — — - S 4—) 6T dz dt 

91 91 C 

0 0 c c 


h t -t h t h 

+ / pcX6T dz + J kX'6T dt - / kX6T' dt 
o oooo 


h t h t 

+ / pc\6T c dz + / S 6a dz. 
0 c O 0 0 


(3.11) 


Note that X(z,t) and s(z,t) are arbitrary functions, and the only 
two unknowns in the above equation are the design derivatives 6T and 
6a. One may now specify the variables X and s in such a way that all 

terms associated with 6T and 6a are dropped. This can be accomplished 
by introducing the following adjoint equations for X and s: 


0 = pc\ + kX~ + pHr |I + S - (2T - 1 J Tdz), 

of oT h 

o 


(3.12) 


o = s + pHrX E + S E 

9a 9a 


(3.13) 


with the terminal conditions. 


X(z,t) = 0, 


0<z<h, 


(3.14) 
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s(z,x) = 0, 0<z<h, (3.15) 

and the boundary conditions, 

X'(0,t) = 0, 0<t<x, (3.16) 

X(h,t) = 0, 0<t<x, (3.17) 

Then, the combination of Eqs. (3.11) - (3.17) provides a simple formula 
for the design derivative of the functional, 

T h . a f af 

64 >= I I (- pcX - pHr — - S — ) 6T c dz dt 

0 0 c c 

h 

+ / pcX(o) 5T (o) dz . (3.18) 

o 

Equation (3.18) shows that the design derivative of , namely, 
6<i>, is a functional of state variables a and T, and the adjoint 
variables X and s. The a and T can be solved by the finite element 
method as discussed in Chapter 2 and Appendix A. The matrix equations 
which solve the nodal values of a and T are mentioned in Eq. (2.14). 
Since the adjoint variables of Eqs. (3.12) - (3.13) form an "adjoint" 
diffusion-reaction system similar to the original one, the same 
numerical scheme used to solve the state variables a and T can be 
extended here to compute the adjoint variables s and X. For instance, 
using the shape functions of a and T in Eqs. (2.9) and (2.11) to 
interpolate the adjoint variables X and s obtains the following matrix 
equations for nodal values of X and s. 

pcCC] X = k[K]X - Q*+ 2 T T [C] - i J T P P T , 


(3.19) 
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CM] S = R* (3.20) 

★ * 

with proper boundary and terminal conditions. The Q and R in the above 
equations have definitions similar to their counterparts of Eqs. (A. 17) 
and (A. 18) defined in Appendix A. 

In general, the adjoint equations can not be solved simultaneously 
with the original system equations. Because of the terminal conditions, 
the adjoint equations can be solved by either the backward integration 
along the real time t - axis directly or the forward integration along 
the artificial time t* - axis, provided that the independent variable t 
is changed to t* as t* = x - t. However, both approaches require the 
solutions of the original system equations prior to solving the adjoint 
equa tions. 


3.3 The Direct Differentiation Method 

The direct differentiation method is an approach that takes 
derivatives of differential equations with respect to a single design 
variable directly. 

For a given functional, Eq. (3.1), and the governing equations of 
the diffusion-reaction system, Eqs. (2.5) and (2.7), the direct 

differentiation results in the following equations in terms of design 

. . .. dT da 

derivatives -rr and — as: 
db db 


x h h h 

<\>' = / [/ 2T~dz--(J Tdz)(J — dz)]dt, 

0 0 


0 


0 


(3.21) 
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dT 


* + 


db 


dT ^ , 

c ,, da 

PC db" + pHr ~adb + pHr 


Sf dT 
dT db 


and 


+ 




da 9f da 3f dT af d c 

db aa db aT db 5T db 

c 


(3.22) 

(3.23) 


It is assumed that T(z,b,t) and a(z,b,t) have enough regularity in 
the time-spatial domain and in the design space. Thus, the order of the 
differentiation is exchangeable, i.e.. 


and 


dT 

db 


dT b 
dt ’ 


da 

db 


da b 
dt ’ 


dT~ 

db 


■ V' 


where the subscript "b" denotes the design derivative. 

Based on the same finite element discretization as used in solving 
the original diffusion-reaction system discussed in Chapter 2, Eqs. 
(3.21) - (3.23) can be converted into a set of matrix equations 

4> b = J T (2T T [C] I b - \ I T P P T I b ) dt, (3.24) 


dT 

pcCC] f = - k[K] T. - pc P — £ + Q, (3.25) 

-b -b - db 

and 

[M] a b = R (3.26) 

The construction of the above equations are discussed in Appendix B. 
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• • 

Note that the coefficient matrices of and in Eqs. (3.25) and 
(3.26) are identical to those of T and a defined in Eqs. (2.14). 
Therefore, the same numerical scheme and numerical tolerance can be 
applied to solve both Eqs. (2.14), (3.25) and (3.26) simultaneously for 
state variables, T and a, and design derivatives, J b and a^. In this 
way, the state variables and the design derivatives achieve the same 
numerical accuracy, though an additional set of equations such as Eqs. 
(3.25) and (3.26) is needed in this approach for each design variable b. 

Regarding the computational efficiency, it is worthwhile mentioning 
two notes here: 

• • 

1. Because the coefficient matrices of T and a are identical to 

-b -b 

• • 

those of T and a, the triangular factorizations of matrices [C] 

and [M] are needed to be done only once. The calculation of 

• • 

T, and a can be carried out by back substitution for each of 
-b -b 

the design variables. 

2. Compared to the original system equations, the right side of 

equations for computing J b and g^, such as Eqs. (3.25) and 
(3.26), may have different frequency contents. Thus, to 
maintain the sa e numerical accuracy, a smaller time 

step At may be required for the DE program to solve the pairs 
(T, a) and (J b 2 b ) simultaneously. 

Once T and T b are available, the design derivative given in Eq. 
(3.24) can be easily obtained by the numerical integration. Another 
suggestion is to rewrite the integral form of Eq. (3.24) as a 
differential equation of given as 
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4^ = 2 I T [C] J b - ^ I T P P T T b . (3.27) 

The above equation of 4> b can then be solved simultaneously with 
equations of (T, ot) and (T^, In this way, one extra design 

derivative <|> b for each design variable is introduced in the design 
sensitivity analysis. However, the accuracy of 4 ^ is secured. Equation 
(3.27) is used in the next section and Chapter 4 to compute the 4 ^. 

3.4 Numerical Examples 

Two examples which deal with the cure process of compression 
molding are presented in this section to discuss the numerical accuracy 
of the adjoint variable method and the direct differentiation method for 
calculating the thermal design derivatives. The accuracy of the thermal 
design sensitivity analysis can be checked by using the fundamental 
definition of design derivatives which can be approximated by the finite 
difference. In other words, it is mathematically true for a small 
perturbation of design variable Ab that 

* = cIF 

= A4> 

“ A &• 

The finite perturbation of the design variable Ab is defined as the 

difference between a perturbed design b* and the nominal design b, i.e., 
★ 

Ab = b - b. As a result, it follows that 

A(|» = <|>(b ) - c|>(b) 


- i>' Ab. 


(3.28) 
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The above equation provides a simple means to check the accuracy of the 
design sensitivity analysis. Nevertheless, the difficulty of this 
method is the selection of Ab. If Ab is too large, the approximation in 
Eq. (3.28) is not valid. On the other hand, if Ab is too small, the 
round-off error in the computation of [ cp( b* ) - 4>(b)] becomes too large 
to apply the approximation of Eq. (3.28) for the sensitivity 

calculation. In the following examples, several values of Ab are used 
to check the accuracy of the design sensitivity analysis in order to 
ensure the quality of the approximation given in Eq. (3.28). 

The first example presented here deals with the cure process Eqs. 
(2.14) in which the cure temperature of the process is assumed to be a 
constant temperature. The nominal cure temperature is taken as 
423°K. According to the approximation defined in Eq. (3.28), the 
results shown in Fig. 3.1 demonstrate that the thermal design 
sensitivity calculated by the direct differentiation method is more 
accurate than the results calculated by the adjoint variable method. 
Moreover, by using the direct differentiation method, one can also get 
the time histories of the design derivatives of state variables as shown 
in Figs. 3.2 and 3.3. 

The second example deals with the same cure process as example 
one. However, the profile of the cure temperature is assumed to be the 
same as the one given in Fig. 3.4 where the heating rate is considered 
as a design variable. The nominal value of the heating rate is taken as 
4.0°K/sec. The resul ts of thermal design sensitivity analyses are shown 
in Fig. 3.5. The time histories of the design derivatives are shown in 
Figs. 3.6 and 3.7. 
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From the above two examples, it is obvious that the direct 
differentiation method provides more accurate results than the adjoint 
variable method does. The direct differentiation method also yields the 
time histories of the design derivatives. In addition, the information 
of design derivatives of the pointwise constraints can be obtained by 
using the direct differentiation method without extra cost. 

Based on the numerical study, it is concluded that the direct 
differentiation method is superior to the adjoint variable technique in 
terms of accuracy and physical interpretation of results. In the next 
chapter, the thermal design sensitivity is to be calculated by using the 
direct differentiation method. 
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CHAPTER IV 

OPTIMAL CURE CYCLE DESIGN 

At the present time, the cure cycle is generally selected by the 
method of parametric study conducted by either physical experiments or 
numerical simulation [1], This "trial and error" type of approach may 
be costly and inefficient. Therefore, the development of a unified and 
systematical approach to design the "best" cure cycle for a given 
composite material is needed. 

The optimal design of cure cycles studied here are to find the 
optimal cure temperature during the cure process to have the 
thermosetting resin cured as uniformly as possible and to assume the 
resin completely cured at the end of the cure process. 

Numerical results presented in this chapter have been obtained by a 

recursive quadratic programming algorithm with an active set strategy, 

called a linearization method [13]. The main reason to use this 

linearization method is that it has been proved to be globally 

convergent. The linearization method, with its active set strategy, 
minimizes the number of constraints which must be considered in each 
design iteration. The detailed derivation of this algorithm can be 
found in reference [13,14]. The concept of this algorithm is that, 
rather than directly solving the optimal criteria, a small perturbation 
for each design variable is determined in each iteration to reduce the 
objective function and correct the violation. Note that, in this 
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approach, the reduction of the objective function and the correction of 

the constraint violations are approximated by the design gradients. In 

general, whenever, the objective function can not be reduced further by 

changing the design variables and all constraints are satisfied, it is 

indicated that the linearization method is converged and the optimal 

solution is obtained. In this algorithm, however, it has been proved 

2 

[15] that a local optimal solution is found when the a norm of the 
perturbation of design variables approaches to zero. The complete flow 
chart of this algorithm is listed in Fig. 4.1. As shown in the figure, 
the proposed CAD method consists of three uncoupled modules: namely, 
optimization, analysis, and sensitivity calculation. Once the analysis 
capability of any cure process is established, the optimization and 
sensitivity calculation modules can be added onto it to constitute a 
unified CAD method and to generate an optimal cure cycle 
systematically. In this thesis a well documented DE program is used for 
the numerical integration of state equations and equations of thermal 
design derivatives. 


4.1 Problem Statement 

As mentioned, optimal design techniques have been successfully 
applied to various transient problems. In general, an optimal design 
problem consists of design variables, the objective, constraint 
functions, and state equations which describe the physical model of 
interest. 

In our study, the process of interest is a compression molding of a 
filled polyester resin reinforced by chopped glass fibers SMC [2]. The 
state equations are limited to a heat conduction equation and an 
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empirical equation which addresses the chemical-kinetic reaction of 
resin. 

Since this study excludes the resin flow and the effects of the 
pressure cycle, the selection of the cure temperature is then limited to 
the following considerations: 

a) The maximum temperature inside the composite during the cure 
process can not exceed a value of 500°K in this study. 

b) At the end of the cure process, the material is cured 
completely. Mathematically, the degree of cure is required to 
be at least 0.85 when the cure process is completed. 

c) The material is cured uniformly at any time during the cure 
process. 

The first two requirements may be formulated as constraints 


: T(t) < 500 °K, 

0<t<T, 

(4.1) 

: <x(t) > 0.85 


(4.2) 


where x is the operational period of the cure process. The x is not 
considered as a design variable in this study. The last requirement can 
be met by formulating it as the objective function of the optimal 
design. Since the degree of cure is a function of temperature, the 
uniform distribution of temperature across the thickness of the material 
implies the cure uniformity. The objective of the optimal design may 
then be set to achieve the temperature uniformity across the section of 
the material during the cure process. The objective function 4> Q is then 
defined as the greatest value of the standard derivation of temperature 
distribution which happens during the cure process 
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4 (T »t) = Maximum [/ T 2 dz - (/ Tdz) 2 /h]. (4.3) 

0<t<x o o 

The minimum of the objective function, <|» q (T ,t) = 0, corresponds to 
a uniform temperature distribution inside the composite at any time 
during the cure process. Consequently, a uniform cure state can be 
accomplished during the complete duration of the cure process. 

To sum up, the optimal cure cycle design is defined as follows: 

"The optimal cure temperature (temperature profile) is designed, 
subject to the limitations regarding the maximum temperature and the 
state of cure, so as to achieve a uniform temperature distribution along 
the cross section of the material at any time during the cure process." 

The mathematical formulation of the stated optimal design problem 
is a mimimax problem given as 

h ? h ? 

Minimize Maximum [/ i dz - (/ Tdz) /h] 

T c (t) 0<t<x o o 

subject to the constraints stated in Eqs. (4.1) - (4.2) where the 
temperature T(t) and the degree of cure a(t) are the solutions of state 
equations. 

It is known that the objective function of a minimax problem is 
discontinuous in the design space [16]. To avoid the computational 
difficulty, one may modify the optimal design problem to a standard form 
by introducing an extra design variable b and an additional constraint 

*3 as 

Minimize = b 
T c (t), b 

subject to 


0<t<x, 


4^ : T(t) < 500 °K, 


(4.4) 
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4»2 • «(x) > 0.85 


(4.5) 


[/ T 2 dz - (/ Tdz) 2 /h] <b, 
o o 


0<t<X. 


(4.6) 


The design function T c (t) can now be parametrized by a linear 
combination of design parameters and given functions such as 

ND * 

T ( t) = Z b. N . (t) 
c 1 ' * 

★ 

where ND is the total number of design variables, b^ , and ( t) are any 
independent functions. In this presentation, the design function is 
simply represented by a combination of linear polynomials and sinusoids 

T (t) = T + b.t + b 0 sin(nt/x) + b-sin(2itt/x) + 

C O i c. 0 


b^sin(3nt/x) , (4.7) 

where T q is the room temperature and bj» b 2 » b 3 , and b^ are to be 
determined by the optimal design algorithm. As a result, the design 
space becomes finite dimensional, and the design sensitivity (or 
gradient) of temperature T(t) and state of cure a(t) are taken as 
derivatives with respect to design parameters bj, b 2 , b 3 , and b^. For 
example, the design derivatives of constraint functions, Eqs. (4.4) - 
(4.6), can be derived as 


"SET 


= ( 


ant) 

ab.. 


)/500. 


(4.8) 


^2 _ / da(x)^ /A oc 

W7 ‘ oFj ) /0,85 » 


(4.9) 



36 


and 


d< N ? h ftT h aT 

ab- 3 b [/ T ab- dz - ( ' Tdz)(/ 5b- dz)/h] - 
i o i o oi 


(4.10) 


It is obvious that the design sensitivities of and are needed 

for the calculation of the above function gradients. Note that the 
pointwise constraints, <|>j and 4^, are imposed at every time grid point 
in this study. The optimal cure cycle T c ( t ) of the above problem can be 
found numerically by calculating the thermal sensitivity analysis and by 
using a gradient-based mathematical programming technique. As discussed 
in the last chapter, the direct differentiation method is employed 
hereafter to compute the required thermal design derivatives. 


4.2 Numerical Examples 

To show the applicability of the proposed CAD scheme for the 
optimal cure cycle design, four examples associated with various problem 
formulations are discussed and presented in this section. 

Example 1. 

The objective function defined by Eq. (4.3) is minimized to find 
an optimal cure cycle design for processing a chopped glass fibers SMC 
with 10 mm in thickness. The total processing time is limited to 100 
seconds. The maximum temperature allowed inside the SMC is 500°K; and 
the degree of cure a is required to reach at least 0.85 at the end of 
the cure process. 

The physical and kinetic properties of the SMC material are given 
in Table 2.1 [2], It requires 22 iterations of the algorithm to obtain 
the optimal solution for this case. Numerical results for the optimal 
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design are listed in Table 4.1. Figure 4.2 shows the cure temperatures 
computed at iterations 1, 11, and 22, respectively. It should be noted 
that the cure cycles designed in this example are always made to start 
from room temperature which is different from the common practice with 
pre-warmed press platens employed in the compression molding process. 
The initial cure temperature (iteration 1) has a peak and a valley at 
around 20 and 60 seconds, respectively. Among the iterations (1, 11, 
and 22) the initial heating rates are seen to decrease gradually; 
eventually a single peak near 70 seconds is obtained for the optimal 
solution. The values of the objective functions (<i> 0 computed during 
iterations) are shown in Fig. 4.3. The optimal cure cycle indeed 
delivers a more uniformly distributed temperature. The optimal solution 
has reduced the maximal value of temperature deviation from 2,700 to 378 
units. 

Temperature distributions across the thickness of the SMC at 
selected instants of 20, 60, and 100 seconds during the cure process are 
compared for iterations 1, 11, and 22, shown in Fig. 4.4a, b, and c, 
respectively. As noted previously in Fig. 4.2, the low initial heating 
rate helps the optimal cure cycle (iteration 22) achieve a more uniform 
temperature distribution at the early 20 seconds into the cure 
process. At the 60 seconds into the cure process, the cure cycle of the 
eleventh iteration obtains a better distribution profile as seen in Fig. 
4.4b. This is also a direct consequence of the negative heating rates 
(cooling) between 30 to 50 seconds observed from Fig. 4.2. On the other 
hand, the high initial heating rate provided by the cure cycle of 
iteration 1 has strongly triggered the kinetic reaction in the first 20 
seconds into the cure process. Consequently, a highly non-uniform 
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temperature is observed between 30 to 50 seconds, even though the cure 
cycle of iteration 1 possesses the steepest cooling rate within the same 
time interval among all of the cure cycles shown in Fig. 4.4. In other 
words, the heat change caused by the steepest cooling rate needs enough 
time to be diffused into the material in order to offset the extra heat 
generated by chemical reactions between 30 to 50 seconds. It is, 
therefore, expected that, as shown in Fig. 4.4c, the more uniform 
temperature distribution can be achieved by the cure cycle of iteration 
1 at the latter stage of the cure process. 

The distributions of the states of cure a across the SMC thickness 
are shown in Figs. 4.5, 4.6, and 4.7 for the cure cycles of iterations 
1, 11, and 22, respectively. In all cases, the cure of the resin starts 
from the surface regime and extends into the core of SMC. The cure 
cycle of iteration 1 requires the shortest time (60 seconds) to complete 
the cure reaction of the material with a > 0.85. 

The low cure temperatures provided by the optimal cure cycle at the 
intial 40 seconds into the cure process can hardly trigger the chemo- 
kinetic reactions of the resin matrix as seen from Fig. 4.7. However, 
for the optimal cure cycle, the rate of the cure process speeds up in 
the next 40 seconds and eventually the degree of cure exceeds the 0.85 
limit throughout the SMC material at the end of the cure cycle. 

Example 2. 

This is identical to the problem stated in Example 1 except that 
the total processing time is extended to 150 seconds. 

The optimal cure cycle of Example 1 is used as the initial trial 
cure cycle studied here. The algorithm needs another 25 iterations to 
reach the optimal solution for this problem. Table 4.2 lists the 
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numerical results for the optimal design. The optimal cure cycles of 

Examples 1 (iteration 1) and 2 (iteration 25) are compared in Fig. 
4.8. It is noted that the cure cycle of Example 2 possesses even slower 
initial heating rates, and the cooling period is no longer needed over 
the entire cure cycle. The values of the objective functions are also 
compared in Fig. 4.9 as indicated by a drop of maximum temperature 
deviation from 250 to 110 units, the extension of the total processing 
time enables an optimally designed cure cycle to achieve better 
temperature uniformity during the cure process. Such improvement is 
also magnified in Fig. 4.10 where temperature distributions across the 
SMC thickness at selected instants during the cure process for the above 
two cure cycles are compared. The distribution of the state of cure for 
the optimal cure cycle of this example is shown in Fig. 4.11. The cure 
behavior basically shows the same characteristics as the optimal 
solution of Example 1. 

Example 3. 

This is identical to the problem stated in Example 1, except that 
the maximum cure state deviations defined by a new objective function 
(Eq. 4.11) is minimized; and the total processing time is extended to 
150 seconds. 

The goal of present cure cycle design is to have the thick 
composite SMC cured uniformly. It is rationalized that perfectly 
uniform temperature distribution across the SMC thickness all the time 
during the cure process will secure such an objective. However, as 
shown by the previous two examples, the optimal cure cycle of Example 2 
does not establish significantly uniform states of cure, even though the 
temperature uniformity (as indicated by the objective functions 4> Q ) has 
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been improved greatly. This could be due to the fact that the degree of 
cure is quite sensitive to the temperature non-uniformity. In this 
example, the objective function Eq. (4.3) and the constraint Eq. (4.6) 
are replaced by the following forms, respectively; 

h ? h ? 

<J> : Maximum [/ oc dz - (/ adz) /h], (4.11) 

o o 

h ? h ? 

(K ; [/ a dz - (/ adz) /h] <b, 0<t<-c. (4.12) 

o o 

The optimal cure cycle of Example 2 is used as the initial trial 
solution in this case. Fifty five more iterations are required to have 
the algorithm converged to the optimal solution. The optimal cure cycle 
shown in Fig. 4.12 has the characteristic of a ramp-hold-ramp feature. 
Numerical results are listed in Table 4.3. The values of the objective 
functions of Eq. (4.11) for the starting and the optimal cure cycles are 
compared in Fig. 4.13. The optimal solution reduces the maximal value 
of the state of cure deviation as expected. However, the differences 
are not at all significant. The distribution of the state of cure, 
shown in Fig. 4.4, corresponding to the optimal cure cycle of Example 3, 
is similar to that of Example 2. The temperature distribution under the 
optimal cure cycle in this example is also shown in Fig. 4.15. These 
results indicate that the change of the objective function from Eq. 
(4.3) to Eq. (4.11) has only minor effects on the uniformity of the 
state of cure, even though the profiles of the optimal cure cycle 
designs are quite different. 
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Example 4. 

This is identical to the problem stated in Example 1, except that 
the maximum state of cure deviation defined by a new objective function 
(Eq. 4.13) is minimized; the total processing time is allowed for 150 
seconds and the initial cure temperature is regarded as an additional 
design variable. 

In this example, the objective function Eq. (4.3) and the 

constraint Eq. (4.6) are replaced by the following forms, respectively; 

<1> : Maximum (a -a ), (4.13) 

0 0<t<t 0 c 

4^: (a Q -a c )^ <b, 0<t<T. (4.14) 

In this example, the objective function becomes the difference of 
the degree of cure between the outer surface, a Q , and the center plane, 
a . The optimal design algorithm does not really achieve 

convergence. However, after the forty-sixth iteration, the changes of 
the design variables become small, compared to the magnitudes of the 
design variables and all the constraints are satisfied. The feasible 
design of the forty-sixth iteration is then regarded as the final 
design. Numerical results of this example are listed in Table 4.4. The 
cure temperature profiles and the values of the objective function of 
Eq. (4.13) for iterations 1, 46, and that of reference [2] are shown in 
Figs. 4.16 and 4.17, respectively. Note that the cure cycle given in 
reference [2] is a constant temperature profile with T c = 423°K. The 
cure temperature profile of the forty-sixth iteration has the 
characteristic of a ramp-hold-ramp feature similar to that of Example 
3. It is recognized that the high onset cure temperature of the optimal 
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cure cycle benefits the cure cycle design by initiating the chemical 
reaction as quick as possible so as to achieve a better state of cure 
distribution. Similar to the previous case, the final solution of this 
example does improve the uniformity of the state of cure, though this 
improvement is insignificant. The final results of the temperature 
distribution and the degree of cure distribution at the selected 
instants are shown in Figs. 4.18 and 4.19, respectively. 
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CHAPTER V 

CONCLUDING REMARKS 

At the present time, the cure cycle is generally selected by the 
method of parametric study conducted by either physical experiments or 
numerical simulation. This type of approach may be costly and 
inefficient. Therefore, the development of a unified and systematical 
approach to design the best cure cycle for a given composite material is 
needed. A unified Computer-Aided Design method is introduced in this 
thesis to design the best cure cycle for a resin-fiber composite 
laminate. The characteristics of curing the resin-fiber composite 

laminate can be simplified as a diffusion-reaction system and a kinetic 
model of the cure rate of the resin. The proposed CAD method consists 
of three uncoupled modules: namely, optimization, analysis, and 

sensitivity calculation. Once the analysis capability of any cure 
process analysis is established, the optimization and sensitivity 
calculation modules can be added onto it to constitute a unified CAD 
method and to generate the optimal cure cycle systematically. 

The optimal cure temperature is required to achieve the following: 

1. The maximum temperature during cure should be less than 500°K. 

2. At the end of the cure process, the cure of the resin should be 
completed. 
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3. The state of cure should be uniform during the cure process. 

It is shown by numerical examples that this method performs satisfac- 
torily by carefully selecting the optimization formulation. 

The finite element discretization is employed to convert the 
initial-boundary value problem of interest into a set of first order 
differential equations. These equations are then solved simultaneously 
by a numerical code called DE. The finite element method provides a 
systematical method to analyze the temperature and the degree of cure 
distribution. Note that very good agreements exist between the finite 
element solutions and the existing solutions. 

Two methods, the adjoint variable technique and the direct 
differentiation method, are studied for the thermal design sensitivity 
analysis. In general, the adjoint equations can not be solved 
simultaneously with the original system equations. Therefore, it is 
difficult for the adjoint variable technique to have the numerical 
accuracy of the design derivatives equal to that of the state 

variables. Numerical examples indicate that the direct differentiation 
method is superior to the adjoint variable technique in terms of 
accuracy and physical interpretation. In this thesis, the direct 

differentiation method is used for the optimization algorithm. 

The optimization algorithm used in this thesis is a recursive 
quadratic programming algorithm with an active set strategy, called a 
linearization method. The main reason to use this linearization method 
is that it has been proved to be globally convergent. The optimal 

design of the cure cycles studied in this thesis is to find the optimal 
cure temperature during the cure process to have the thermosetting resin 
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cured as uniformly as possible and to secure the cure completeness of 
the resin at the end of the cure process. 

Four optimal design formulations for the optimal cure cycle design 
of chopped-fiber SMC are studied. The first example deals with the 
minimum temperature deviation across the thickness during the cure 
process, and the total process time is 100 seconds. The optimal cure 
cycle of this example reduces the value of the objective function from 
the initial design 2,700 to the optimal one 378 units. The second 
example deals with the same conditions as Example 1 except that the 
total processing time is extended to 150 seconds. Using the optimal 
solution of Example 1 as the initial design, the maximum temperature 
deviation drops from 250 to 110 units. The third example is identical 
to Example 1, except that the objective function is changed to minimize 
the degree of cure deviation across the thickness along the cure process 
and that the total processing time is 150 seconds. The optimal solution 
reduces the value of the objective function from the initial design 
0.04691 to the optimal one 0.03385. The last example deals with the 
minimization of the difference between the degree of cure of the outer 
surface and the center plane. The total processing time is 150 seconds 
and the initial cure temperature is regarded as an additional design 
variable. This example is different from the previous examples in which 
the initial temperatures are set to be the room temperature. The 
optimal cure cycle reduces the degree of cure deviation between the 
outer surface and the center plane from the initial design 0.8032 to a 
feasible design 0.6992, and increases the initial cure temperature to 
395° K. Comparing the optimal cure cycle with that of reference [2], 
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the optimal cure cycle does make the state of cure distributed more 
uniformly during the cure process. 

The numerical examples show that the proposed CAD approach is 
valid. It does improve the degree of cure uniformity. This study also 
indicates that the total processing time, the initial heating rate and 
the initial cure temperature have significant effects on the cure 
uniformity. Enlarging the design space by increasing the total 
processing time and adding more design variables can further improve the 
cure cycle design. The long CPU time is a concern in the optimization 
process. To shorten the CPU time one may use other numerical 
algorithms, such as an unconditional stable algorithm, which can 
increase the time step, or the other kind of optimization algorithms, 
such as the I-DESIGN [17], which does not need line search to accelerate 
the optimization process. 

The proposed CAD scheme can be extended to design an optimal cure 
cycle for the autoclave processing. Note that, to design a realistic 
cure process of resin-fiber composite optimally, some modification over 
the presented simplified SMC model are required. For example, the flow 
model, the chemo-viscosi ty and the heat convection should be considered 
in the analytical model. The total processing time should be regarded 
as a design variable to enlarge the design space. Moreover, some 
experiments are necessary to verify the results of the cure cycle 
designs. 

In summary, it is clear and convincing that the optimization is a 
very promising technique to enhance the capability of designing the 
optimal cure cycle for the cure process of resin-fiber composite 
laminate. 
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APPENDIX A 


FINITE ELEMENT MATRIX EQUATIONS 
FOR A DIFFUSION-REACTION SYSTEM 


The equations (2.12) - (2.13) can be rewritten here: 

NE Z i+1 

0=2 / (pcTu + kT'u' + pcT u - pHrfu) dz, (A.l) 

1 z. C 

l 

and 

NE z i+l 

0 = E / (ct-f) v dz (A. 2) 

1 z i 

where NE denotes the number of elements and u(z) and v(z) are the 
testing functions. 

According to the Galerkin method, the arbitrary interpolation 
functions Nj and Lj defined in Eqs. (-2.8) and (2.10) can be considered 
as u(z) and v(z), respectively. Furthermore, the temperature and the 
degree of cure can be expressed as linear combinations of interpolation 
functions and the nodal values as shown in Eqs. (2.9) and (2.11). One 
may then write a matrix equation for the first integral of Eq. (A.l) 
defined as 
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. i+l • 

/ pcT u dz 

i 

X. 

1 -r 

= / pc N. N dz T, j=l,2,3,4. (A. 3) 

o J 


where z denotes the global coordinate system, z denotes the local 
coordinate system, z - z-z.. , and JL = z^-z... Let a vector C be 
defined as 

X . 

T 1 t - 

C = / N. N dz, j = 1,2, 3, 4. 

’j o J " 


Consequently, one has 


/ pc T u dz 
z . 


= P c C T T (A. 4) 

j 


for j = 1 to 4. At this moment, one may introduce an element 
matrix [C] 1 in the following way: 


CC] 1 = / 


*i 


N N T dz 


= [C , C , C , C ] T 
12 3 4 

where the vector C is defined in Eq. (A. 4). It is understood that 
1 ~ j 

[C] is symmetric. Thus, for u = Nj with j = 1 to 4, one has 


= pc CC] T 


(A. 5) 


Similarly, using the following definitions for an element matrix 
• • 

[K] 1 and a vector P 1 : 


and 


[K] 1 = / r r T dz, 
o 

X . 

.• 1 
P = / N dz 
o 


(A. 6) 


(A. 7 ) 


one obtains matrix equations for the second and third integrals of Eq. 
(A.l) as 

z i + l 

/ k T' u' dz 

z . 

1 


and 


= k [K] 1 ' I 



/ pc T u dz 
c 

z. 


= pc 



(A. 8) 


(A. 9) 


The element matrices [K] 1 , [C] 1 and P 1 can be expressed explicitly in 
terms of the length of element i, X. , as follows: 
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[K] 1 = 


51 . 

i 


1_ 

10 


21 . 

TF~ 


Sym. 


6 

5*. 


1 

Iff 


6 

HT 


-1 

10 


*i 


1 

1U 

21. 

1 

IT 


(A. 10) 


[C] 1 = 


and 


l. 

l 

TZU 


156 22 A. 54 13JL 


4JL^ 13A. -3A? 


Sym. 156 -22A ‘ 




P’ = 


2 2 
{*./2, A./12, A./2, -A^/12} 


(A. 11) 


(A. 12) 


Following the previous derivation, for v = Lj with j = 1, 2, one 
may write a matrix equation for the first integral of Eq. {A. 2) as 


z i+l 

/ av dz 
z. 


= / L, L T dz a 

n J ~ 


(A. 13) 
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where Lj is a shape function for the degree of cure, 
matrix [M] 1 be defined as 

[M]' - J L L T dz. 

0 


Let an element 


(A. 14) 


Thus, for v = Lj with j=l , 2, one has the following identity: 


Z i+1 


/ a v dz 
z i 

= [M] 1 a. 


(A. 15) 


The matrix [M] for the ith element can be derived explicitly in terms of 


JL as 


CM ] 1 = 


Aj/3 A./6 


A 1 /6 A./3 


(A. 16) 


After substituting an interpolation function for the trial function 
u(z), the last term in Eq. (A.l) yields a scalar Q^p for ith element as 


A i 

Q ip = / pHrf (a,T,T c ) N p (z) dz. 


p = 1,2, 3, 4. (A. 17) 


In a similar fashion, the discretized form of the last term in Eq. (A. 2) 
can be expressed as 



*i 

I f(a,T,T ) L (z) dz, 
o C q 


q = 1,2. 


(A. 18) 
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The cure rate a = f(a,T,Tc) in Qjp and R-jg is a function of the 
temperature and the degree of cure defined in Eq. (2.5). Note that the 
temperature distribution has been replaced by T+T c . Thus, Eq. (2.5) 
becomes 

a = f(a,T,T c ) 


= Uje 


J 1/R(T+T C ) 


-d 


a^e 


2/R(T+T ) 

C a m )(l-ct) n . (A. 19) 


Within the ith element, i.e., z. < z < z^, the distributions of 
temperature is interpolated by shape functions and nodal variables as 
defined in Eq. (2.9). To abbreviate the notations, the following r is 
introduced: 

r = T(z,t) + T c (t) 


= N T I + T c (t) (A. 20) 

It is also understood that 

a(z,t) = L T a. (A. 21) 

Hence, Q^p and R^g can be expressed as 


Q. = pHr f 1 U 1 e dl/RF + a 2 e d2/Rr <x m )(l-a) n N (z)dz, p = 1,2, 3, 4. 
• o ^ 


and 


o r V / d i/Rr x a / d 2/Rr m wi >n . , , 

R. = J (a^e + a^e a )(l-a) L (z)dz. 


q = 1.2 


where Q^p and are components of a 4 x 1 load vector Q and a 2 x 1 
load vector R, respectively. Note that it is difficult to evaluate the 
above equations analytically because of the complexity. However, using 
Eqs. (A. 20) - (A. 21), the Q^p and R^q can be computed numerically 

provided that the nodal values of T and a are known. The numerical 

examples represented in Section 2.3 are solved by using the Simpson's 
rule for this numerical integration. The element matrix equations of 
Eqs. (A.l) and (A. 2) for the ith element can then be shown as 


pcCC] 1 T + kCK] 1 T + pc P 1 T c = Q 1 , 


(A. 22) 


[M : 1 a = R 1 , 


(A. 23) 


Finally, these element matrix equations can be assembled into a set 
of global system equations in a matrix form as 


pc[C] T + k[K] T + pc P T = Q 



APPENDIX B 


FINITE ELEMENT MATRIX EQUATIONS 
OF THERMAL DESIGN SENSITIVITY ANALYSIS 
BY USING THE DIRECT DIFFERENTIATION METHOD 


The given equation (A. 19) can be rewritten here: 
« = f(«.T.T c ) 


-d l/R(T+T ) ~ d 2/R(T+T ) 

■ (a x e + a 2 e c a m )(l-a) n . (B.l) 

Following the definition of the direct differentiation method, one 
takes derivatives of the above differential equation with respect to a 
single design variable directly to obtain 

da df 

m = m 


_ 3f da 5f dT 5f dT c 

da db + 9T db + 5T db“* 

c 

Since f is an explicit function of a, T, and T c , one can find 
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5f 

5a 


-d 


-n(a x ( 


1/R(T+T ) 


-d 


2/R(T+T ) 


+ a 2 e 


a m )(l-a) 


n-1 


and 


-d 2/R(T+T ) 


m-1 


a 


(l-a) n 


(B .3 ) 


5f 

Fr 


R(T+T ) 

C 


’ d l/R(T+T ) “ d 2/R(T+T ) 

9 (a^d^e + a 2 d 2 6 C am )d“ a ) n 


5f 

Fr 


the equations (3.22) and 3.23) can be rewritten here as 


PC f b 


and 


a 


b 


kT~ - pcf + pHr 4“ a k + P Hr T + P Hr 
b cb 5a b 5T b 


T 

5T cb’ 


5f 5f _ 5f _ 

— a + — T + — T 
9a b 9Tb 9T cb 

c 


(B.5) 


(B.6) 


It is assumed that T b (z,b,t) and a b (z,b,t) have the same regularity 
as T(z,b,t) and a(z,b,t). Thus, the shape functions of T(z,b,t) and 

a(z,b,t) can be employed here to interpolate the T b (z,b,t) and 

a(z,b,t), i.e. 


T b (z,b,t) = N T (z) T b , 
and 

a b (z,b,t) = L T (z) a^. 


(B.7) 

(B.8) 


where N(z) and L(z) are defined in Eqs. (2.8) and (2.10), and J b and 
2 b are the vectors of nodal values of T b (z,b,t) and a b (z,b,t). 
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Based on the same derivation discussed in Appendix A, the global 
matrix equations for Eqs. (B.5) and (B.6) can be shown as 

pc [C] T + k[K] T + pc P T = Q, ( B .9 ) 

-d -b cb 

and 

M] a b = R (B.10) 


where matrices [C], [K], [M], and P are the same as those given in Eqs. 
( A.10)-(A.13) and (A. 16), and the components of Q and R for the ith 
element are evaluated as 



l. 

i 


= pHr / '( 


0 


a f af T af 

— a + — T + 

da b 5T b 5T 

c 


T JN (z)dz, 
cb' p 


and 



f 5f 5f 

l 5a b 5T b 



T JL (z)dz, 
cb y q 


p = 1,2, 3, 4. (B.ll) 


q = 1,2. (B.12) 


The N(z) and L(z) in the above equations are the shape functions which 
interpolate the design derivatives of the temperature and the state of 
cure, respectively. Note that it is difficult to calculate Eqs. (B.ll)- 
(B.12) analytically because of complexity. However, the and 
can be computed numerically provided that the nodal values of T, a, T , 
and are known. 
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Table 4.1 Convergence History of Optimum 
Design (Example 1) 



Initial Value 

Final Value 

Design Variables: 

b l 

100. 

115.74 

b 2 

100. 

100.52 

b 3 

100. 

-5.16 

b 4 

100. 

6.83 

cost function 

2,700. 

378.33 

2 

H norm 

33.62 

0.0054 

Computer Time 

_ _ __ 

108,140 CPU sec. 

(on CDC 855) 

No. of Iterations 

— 

22 
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Table 4.2 Convergence History of Optimum 
Design (Example 2) 



Initial Value 

Final Value 

Design Variables: 

b l 

115.74 

162.32 

b 2 

100.52 

30.98 

b 3 

-5.16 

-0.18 

b 4 

6.83 

5.00 

cost function 

250.00 

110.86 

2 

l norm 

40.88 

0.0147 

Computer Time 

* V — 

156.370 CPU sec. 

(on CDC 855) 

No. of Iterations 

— 

25 
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Table 4.3 Convergence History of Optimum 
Design (Example 3) 


Initial Value Final Value 


Design Variables: 


b l 

162.32 

200.42 

b 2 

30.98 

25.04 

b 3 

-0.18 

46.99 

b 4 

5.00 

-5.28 

cost function x 10^ 

46.91 

33.85 

2 

Z norm 

5.29 

0.19 

Computer Time 
(on CDC 855) 

— 

319,000 CPU sec 

No. of Iterations 

— 
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Design Variables: 
To 


b l 

b 2 

b 3 

b 4 

cost function 
2 

Z norm 

Computer Time 
(on IBM 4381) 


No. of Iterations 


300.00 


395.85 


168.39 

106.59 

47.41 

-46.92 

40.67 

23.42 

-0.99 

-24.29 

0.8032 

0.6992 

111.40 

4.08788 



196,800 CPU sec. 
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Finite Clement Results 



Figure 2.2 One Dimensional Heat Conduction with Prescribed 

End Temperature: Spatial Temperature Distribution. 
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Figure 2.4 Cure Profiles for 10-mm Thick Sheet. 
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Figure 2.5 Temperature Profile at the Center of 
Laminate for Thick Section Autoclave 
Cured Composite. 
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Direct Differentiation 


Adjoint Variable 



Figure 3.1 Thermal Design Derivatives for Press Molding 
with Respect to the Mold Temperature. 


Thermal Design Derivative dT/dT 
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Figure 3.2 Profiles of Thermal Design Derivatives of 

Temperature with Respect to the Mold Temperature. 




Thermal Design Derivative dOr/dTc 
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Direct Differentiation 
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Figure 3.3 Profiles of Thermal Design Derivatives of the State 
of Cure with Respect to the Mold Temperature. 



Change of Functional, 


76 


Direct Different! a t i o n 




Adjoint Variable. 



Change of Design Variable 


Figure 3.5 Thermal Design Derivatives for Press Molding with 
Respect to the Heating Rate. 


Thermal Design Derivative dT/dT 
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Direct Differentiation 


Finite Difference 
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Figure 3.6 Profiles of Thermal Design Derivatives of 

Temperature with Respect to the Heating Rate 
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Figure 4.1 Flow Chart of the Proposed CAD Method. 
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Figure 4.2 Example 1: Profiles of the Cure Temperature 



Iteration 11 



Figure 4.3 Example 1: Profiles of the Least-Squared Temperature Deviation 



Iteration II 



Figure 4.4 Example 1: Profiles of Temperature Distribution. 
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Figure 4.7 Example 1: The Degree of Cure Distribution at 
Iteration 22. 
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Figure 4.8 Example 2: Profiles of the Cure Temperature 
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Iteration 25 



Figure 4.10 Example 2: Profiles of the Temperature Distribution. 
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Figure 4.13 Example 3: Profiles of the Least-Squared Degree of Cure Deviation 
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Figure 4.16 Example 4: Profiles of the Cure Temperature 
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Example 4: The Temperature Distr 
Iteration 46. 






